This notebook runs a longitudinal-functional analysis on simulated data.

library(tidyverse)
library(MASS)
library(Rcpp)
library(splines)
library(LFBayes)
library(plotly)
### Control parameters

set.seed(999)
source("/Users/johnshamshoian/Documents/R_projects/LFBayes/Rfuns/Example/Simulation_funcs.R")
errorvar <- .025
SS <- 20
TT <- 20
t <- seq(from = 0, to = 1, length.out = TT)
s <- seq(from = 0, to = 1, length.out = SS)
n <- 60
tt <- list()
tt[[1]] <- 1:(TT*SS)
tt <- rep(tt, n)
p1 <- 12
p2 <- 12
q1 <- 3
q2 <- 3
Bt <- bs(t, df = p1, intercept = TRUE)
Bs <- bs(s, df = p2, intercept = TRUE)
Bt1 <- bs(t, df = p1, intercept = TRUE)
Bs1 <- bs(s, df = p2, intercept = TRUE)
### Generate key model quantities

H <- GenerateH(q1, q2)
mu1 <- GenerateMu1(s,t)
Lambda <- Loading.Matern(t, p1, q1, Bt)
Gamma <- Loading.Brown.Bridge(s, p2, q2)
Cov <- kronecker(Bs%*%Gamma, Bt%*%Lambda)%*%H%*%t(kronecker(Bs%*%Gamma, Bt%*%Lambda)) + errorvar * diag(SS * TT)
### Generate data from model

x <- mvrnorm(n, mu  = as.vector(mu1), Sigma = Cov)
sx <- sd(x)
mx <- mean(x)
x <- (x-mx)/sx
Smooth_scaled_cov <- (Cov - errorvar * diag(SS * TT)) / sx^2
mu <- (mu1 - mx)/sx
y <- lapply(1:n, function(i) x[i,])
missing <- list()
for(ii in 1:n){
  missing[[ii]] <- numeric(0)
}
X <- cbind(rep(1, n))
### Visualize a few trajectories

nsub <- 6
small_data <- tibble(id = numeric(),
                     func_time = numeric(),
                     long_time = numeric(),
                     value = numeric())
small_data <- small_data %>% 
  add_row(id = rep(1:nsub, each = SS * TT),
          func_time = rep(rep(t, SS), nsub),
          long_time = rep(rep(s, each = TT), nsub),
          value = c(t(x[1:nsub,])))

id.labs <- paste("Subject", 1:nsub)
names(id.labs) <- "1":nsub
small_data %>%
  filter(long_time %in% c(s[5], s[10], s[15])) %>%
  ggplot(aes(func_time, value)) +
  geom_line(aes(linetype = factor(long_time))) +
  facet_wrap(. ~ id, labeller = labeller(id = id.labs)) + 
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), 
        strip.background = element_rect(
        color="black", fill="#FFFFFF")) + 
  labs(x = "Functional time", y = "Response", linetype = "Longitudinal time") +
  scale_linetype_discrete(labels = round(c(s[5], s[10], s[15]), 2))

### MCMC control parameters

iter <- 100 # Number of iterations
burnin <- 10 # Burnin iterations
thin <- 1 # Thinning for each chain
nchain <- 1 # Number of chains
q1s <- 3 # Number of latent factors for functional dimension
q2s <- 3 # Number of latent factors for longitudinal dimension
alpha <- .05 # Type 1 error
neig <- 3 # Number of eigenfunctions for inference
### Processing

MCMC <- mcmcWeakChains(y, missing, X, Bs1, Bt1, q1s, q2s, iter, thin, burnin, nchain, 1)
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
MCMC_eigen <- eigenLFChains(Bs1, Bt1, MCMC, neig, iter, burnin, nchain, s, t, alpha)
### Prepare data for plotting
postmean <- matrix(MCMC_eigen$postmean, nrow = TT, ncol = SS)
fig <- plot_ly(z = ~postmean, x = s, y = t) %>% add_surface()
eigenfunction_tibble <- tibble(time = numeric(), type = character(),
                               value = numeric(), bound = character(),
                               number = numeric())
eigenfunction_tibble <- eigenfunction_tibble %>%
  add_row(value = c(MCMC_eigen$eigvecFuncmean,
                    MCMC_eigen$eigvecFunclower,
                    MCMC_eigen$eigvecFuncupper),
          number = rep(rep(1:neig, each = TT), 3),
          bound = rep(c("mean", "lower", "upper"), each = neig * TT),
          type = "Functional",
          time = rep(t, 3 * neig)) %>%
  add_row(value = c(MCMC_eigen$eigvecLongmean,
                    MCMC_eigen$eigvecLonglower,
                    MCMC_eigen$eigvecLongupper),
          number = rep(rep(1:neig, each = SS), 3),
          bound = rep(c("mean", "lower", "upper"), each = neig * SS),
          type = "Longitudinal",
          time = rep(s, 3 * neig))
### Plot mean surface

axx <- list(
  title = "Longitudinal time"
)

axy <- list(
  title = "Functional time"
)

axz <- list(
  title = "Posterior mean"
)

fig <- plot_ly() %>%
  add_surface(z = ~postmean, x = s, y = t) %>%
  layout(scene = list(xaxis=axx,yaxis=axy,zaxis=axz, aspectratio = list(x = .9, y = .8, z = 1))) %>%
  hide_colorbar()
fig

NA
### Plot eigenfunctions

options(repr.plot.width=6, repr.plot.height=3)
number.labs <- paste("Eigenfunction", 1:neig)
names(number.labs) <- c("1":neig)
eigenfunction_tibble %>%
  ggplot(aes(time, value)) +
  geom_line(aes(linetype = bound)) +
  facet_grid(type ~ number, labeller = labeller(number = number.labs)) +
  scale_linetype_manual(values=c("dashed", "solid", "dashed"))+
  theme_bw() +
  theme(legend.position = "none", strip.background = element_rect(
        color="black", fill="#FFFFFF"), 
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

NA
### Percent variability explained by first few eigenfunctions in functional direction
100 * rowMeans(MCMC_eigen$eigvalFunc)[p1:(p1-neig+1)]
[1] 56.731567 27.272872  7.160911
### Percent variability explained by first few eigenfunctions in longitudinal direction
100 * rowMeans(MCMC_eigen$eigvalLong)[p2:(p2-neig+1)]
[1] 71.550694 16.300779  4.075942
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICBkZl9wcmludDogcGFnZWQKLS0tCgpUaGlzIG5vdGVib29rIHJ1bnMgYSBsb25naXR1ZGluYWwtZnVuY3Rpb25hbCBhbmFseXNpcyBvbiBzaW11bGF0ZWQgZGF0YS4KCmBgYHtyLCBtZXNzYWdlID0gRkFMU0V9CmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KE1BU1MpCmxpYnJhcnkoUmNwcCkKbGlicmFyeShzcGxpbmVzKQpsaWJyYXJ5KExGQmF5ZXMpCmxpYnJhcnkocGxvdGx5KQpgYGAKCmBgYHtyfQojIyMgQ29udHJvbCBwYXJhbWV0ZXJzCgpzZXQuc2VlZCg5OTkpCnNvdXJjZSgiL1VzZXJzL2pvaG5zaGFtc2hvaWFuL0RvY3VtZW50cy9SX3Byb2plY3RzL0xGQmF5ZXMvUmZ1bnMvRXhhbXBsZS9TaW11bGF0aW9uX2Z1bmNzLlIiKQplcnJvcnZhciA8LSAuMDI1ClNTIDwtIDIwClRUIDwtIDIwCnQgPC0gc2VxKGZyb20gPSAwLCB0byA9IDEsIGxlbmd0aC5vdXQgPSBUVCkKcyA8LSBzZXEoZnJvbSA9IDAsIHRvID0gMSwgbGVuZ3RoLm91dCA9IFNTKQpuIDwtIDYwCnR0IDwtIGxpc3QoKQp0dFtbMV1dIDwtIDE6KFRUKlNTKQp0dCA8LSByZXAodHQsIG4pCnAxIDwtIDEyCnAyIDwtIDEyCnExIDwtIDMKcTIgPC0gMwpCdCA8LSBicyh0LCBkZiA9IHAxLCBpbnRlcmNlcHQgPSBUUlVFKQpCcyA8LSBicyhzLCBkZiA9IHAyLCBpbnRlcmNlcHQgPSBUUlVFKQpCdDEgPC0gYnModCwgZGYgPSBwMSwgaW50ZXJjZXB0ID0gVFJVRSkKQnMxIDwtIGJzKHMsIGRmID0gcDIsIGludGVyY2VwdCA9IFRSVUUpCmBgYAoKYGBge3J9CiMjIyBHZW5lcmF0ZSBrZXkgbW9kZWwgcXVhbnRpdGllcwoKSCA8LSBHZW5lcmF0ZUgocTEsIHEyKQptdTEgPC0gR2VuZXJhdGVNdTEocyx0KQpMYW1iZGEgPC0gTG9hZGluZy5NYXRlcm4odCwgcDEsIHExLCBCdCkKR2FtbWEgPC0gTG9hZGluZy5Ccm93bi5CcmlkZ2UocywgcDIsIHEyKQpDb3YgPC0ga3JvbmVja2VyKEJzJSolR2FtbWEsIEJ0JSolTGFtYmRhKSUqJUglKiUKICB0KGtyb25lY2tlcihCcyUqJUdhbW1hLCBCdCUqJUxhbWJkYSkpICsgZXJyb3J2YXIgKiBkaWFnKFNTICogVFQpCmBgYAoKYGBge3J9CiMjIyBHZW5lcmF0ZSBkYXRhIGZyb20gbW9kZWwKCnggPC0gbXZybm9ybShuLCBtdSAgPSBhcy52ZWN0b3IobXUxKSwgU2lnbWEgPSBDb3YpCnN4IDwtIHNkKHgpCm14IDwtIG1lYW4oeCkKeCA8LSAoeC1teCkvc3gKU21vb3RoX3NjYWxlZF9jb3YgPC0gKENvdiAtIGVycm9ydmFyICogZGlhZyhTUyAqIFRUKSkgLyBzeF4yCm11IDwtIChtdTEgLSBteCkvc3gKeSA8LSBsYXBwbHkoMTpuLCBmdW5jdGlvbihpKSB4W2ksXSkKbWlzc2luZyA8LSBsaXN0KCkKZm9yKGlpIGluIDE6bil7CiAgbWlzc2luZ1tbaWldXSA8LSBudW1lcmljKDApCn0KWCA8LSBjYmluZChyZXAoMSwgbikpCmBgYAoKYGBge3IsIGZpZy53aWR0aD04LCBmaWcuaGVpZ2h0PTR9CiMjIyBWaXN1YWxpemUgYSBmZXcgdHJhamVjdG9yaWVzCgpuc3ViIDwtIDYKc21hbGxfZGF0YSA8LSB0aWJibGUoaWQgPSBudW1lcmljKCksCiAgICAgICAgICAgICAgICAgICAgIGZ1bmNfdGltZSA9IG51bWVyaWMoKSwKICAgICAgICAgICAgICAgICAgICAgbG9uZ190aW1lID0gbnVtZXJpYygpLAogICAgICAgICAgICAgICAgICAgICB2YWx1ZSA9IG51bWVyaWMoKSkKc21hbGxfZGF0YSA8LSBzbWFsbF9kYXRhICU+JSAKICBhZGRfcm93KGlkID0gcmVwKDE6bnN1YiwgZWFjaCA9IFNTICogVFQpLAogICAgICAgICAgZnVuY190aW1lID0gcmVwKHJlcCh0LCBTUyksIG5zdWIpLAogICAgICAgICAgbG9uZ190aW1lID0gcmVwKHJlcChzLCBlYWNoID0gVFQpLCBuc3ViKSwKICAgICAgICAgIHZhbHVlID0gYyh0KHhbMTpuc3ViLF0pKSkKCmlkLmxhYnMgPC0gcGFzdGUoIlN1YmplY3QiLCAxOm5zdWIpCm5hbWVzKGlkLmxhYnMpIDwtICIxIjpuc3ViCnNtYWxsX2RhdGEgJT4lCiAgZmlsdGVyKGxvbmdfdGltZSAlaW4lIGMoc1s1XSwgc1sxMF0sIHNbMTVdKSkgJT4lCiAgZ2dwbG90KGFlcyhmdW5jX3RpbWUsIHZhbHVlKSkgKwogIGdlb21fbGluZShhZXMobGluZXR5cGUgPSBmYWN0b3IobG9uZ190aW1lKSkpICsKICBmYWNldF93cmFwKC4gfiBpZCwgbGFiZWxsZXIgPSBsYWJlbGxlcihpZCA9IGlkLmxhYnMpKSArIAogIHRoZW1lX2J3KCkgKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gOTAsIHZqdXN0ID0gMC41LCBoanVzdD0xKSwgCiAgICAgICAgc3RyaXAuYmFja2dyb3VuZCA9IGVsZW1lbnRfcmVjdCgKICAgICAgICBjb2xvcj0iYmxhY2siLCBmaWxsPSIjRkZGRkZGIikpICsgCiAgbGFicyh4ID0gIkZ1bmN0aW9uYWwgdGltZSIsIHkgPSAiUmVzcG9uc2UiLCBsaW5ldHlwZSA9ICJMb25naXR1ZGluYWwgdGltZSIpICsKICBzY2FsZV9saW5ldHlwZV9kaXNjcmV0ZShsYWJlbHMgPSByb3VuZChjKHNbNV0sIHNbMTBdLCBzWzE1XSksIDIpKQpgYGAKYGBge3J9CiMjIyBNQ01DIGNvbnRyb2wgcGFyYW1ldGVycwoKaXRlciA8LSA1MDAwICMgTnVtYmVyIG9mIGl0ZXJhdGlvbnMKYnVybmluIDwtIDEwMDAgIyBCdXJuaW4gaXRlcmF0aW9ucwp0aGluIDwtIDEgIyBUaGlubmluZyBmb3IgZWFjaCBjaGFpbgpuY2hhaW4gPC0gMSAjIE51bWJlciBvZiBjaGFpbnMKcTFzIDwtIDMgIyBOdW1iZXIgb2YgbGF0ZW50IGZhY3RvcnMgZm9yIGZ1bmN0aW9uYWwgZGltZW5zaW9uCnEycyA8LSAzICMgTnVtYmVyIG9mIGxhdGVudCBmYWN0b3JzIGZvciBsb25naXR1ZGluYWwgZGltZW5zaW9uCmFscGhhIDwtIC4wNSAjIFR5cGUgMSBlcnJvcgpuZWlnIDwtIDMgIyBOdW1iZXIgb2YgZWlnZW5mdW5jdGlvbnMgZm9yIGluZmVyZW5jZQpgYGAKCmBgYHtyfQojIyMgUHJvY2Vzc2luZwoKTUNNQyA8LSBtY21jV2Vha0NoYWlucyh5LCBtaXNzaW5nLCBYLCBCczEsIEJ0MSwKICAgICAgICAgICAgICAgICAgICAgICBxMXMsIHEycywgaXRlciwgdGhpbiwgYnVybmluLCBuY2hhaW4pCk1DTUNfZWlnZW4gPC0gZWlnZW5MRkNoYWlucyhCczEsIEJ0MSwgTUNNQywgbmVpZywKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGl0ZXIsIGJ1cm5pbiwgbmNoYWluLCBzLCB0LCBhbHBoYSkKYGBgCgpgYGB7cn0KIyMjIFByZXBhcmUgZGF0YSBmb3IgcGxvdHRpbmcKcG9zdG1lYW4gPC0gbWF0cml4KE1DTUNfZWlnZW4kcG9zdG1lYW4sIG5yb3cgPSBUVCwgbmNvbCA9IFNTKQpmaWcgPC0gcGxvdF9seSh6ID0gfnBvc3RtZWFuLCB4ID0gcywgeSA9IHQpICU+JSBhZGRfc3VyZmFjZSgpCmVpZ2VuZnVuY3Rpb25fdGliYmxlIDwtIHRpYmJsZSh0aW1lID0gbnVtZXJpYygpLCB0eXBlID0gY2hhcmFjdGVyKCksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB2YWx1ZSA9IG51bWVyaWMoKSwgYm91bmQgPSBjaGFyYWN0ZXIoKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG51bWJlciA9IG51bWVyaWMoKSkKZWlnZW5mdW5jdGlvbl90aWJibGUgPC0gZWlnZW5mdW5jdGlvbl90aWJibGUgJT4lCiAgYWRkX3Jvdyh2YWx1ZSA9IGMoTUNNQ19laWdlbiRlaWd2ZWNGdW5jbWVhbiwKICAgICAgICAgICAgICAgICAgICBNQ01DX2VpZ2VuJGVpZ3ZlY0Z1bmNsb3dlciwKICAgICAgICAgICAgICAgICAgICBNQ01DX2VpZ2VuJGVpZ3ZlY0Z1bmN1cHBlciksCiAgICAgICAgICBudW1iZXIgPSByZXAocmVwKDE6bmVpZywgZWFjaCA9IFRUKSwgMyksCiAgICAgICAgICBib3VuZCA9IHJlcChjKCJtZWFuIiwgImxvd2VyIiwgInVwcGVyIiksIGVhY2ggPSBuZWlnICogVFQpLAogICAgICAgICAgdHlwZSA9ICJGdW5jdGlvbmFsIiwKICAgICAgICAgIHRpbWUgPSByZXAodCwgMyAqIG5laWcpKSAlPiUKICBhZGRfcm93KHZhbHVlID0gYyhNQ01DX2VpZ2VuJGVpZ3ZlY0xvbmdtZWFuLAogICAgICAgICAgICAgICAgICAgIE1DTUNfZWlnZW4kZWlndmVjTG9uZ2xvd2VyLAogICAgICAgICAgICAgICAgICAgIE1DTUNfZWlnZW4kZWlndmVjTG9uZ3VwcGVyKSwKICAgICAgICAgIG51bWJlciA9IHJlcChyZXAoMTpuZWlnLCBlYWNoID0gU1MpLCAzKSwKICAgICAgICAgIGJvdW5kID0gcmVwKGMoIm1lYW4iLCAibG93ZXIiLCAidXBwZXIiKSwgZWFjaCA9IG5laWcgKiBTUyksCiAgICAgICAgICB0eXBlID0gIkxvbmdpdHVkaW5hbCIsCiAgICAgICAgICB0aW1lID0gcmVwKHMsIDMgKiBuZWlnKSkKCgpgYGAKCmBgYHtyfQojIyMgUGxvdCBtZWFuIHN1cmZhY2UKCmF4eCA8LSBsaXN0KAogIHRpdGxlID0gIkxvbmdpdHVkaW5hbCB0aW1lIgopCgpheHkgPC0gbGlzdCgKICB0aXRsZSA9ICJGdW5jdGlvbmFsIHRpbWUiCikKCmF4eiA8LSBsaXN0KAogIHRpdGxlID0gIlBvc3RlcmlvciBtZWFuIgopCgpmaWcgPC0gcGxvdF9seSgpICU+JQogIGFkZF9zdXJmYWNlKHogPSB+cG9zdG1lYW4sIHggPSBzLCB5ID0gdCkgJT4lCiAgbGF5b3V0KHNjZW5lID0gbGlzdCh4YXhpcz1heHgseWF4aXM9YXh5LHpheGlzPWF4eiwgYXNwZWN0cmF0aW8gPSBsaXN0KHggPSAuOSwgeSA9IC44LCB6ID0gMSkpKSAlPiUKICBoaWRlX2NvbG9yYmFyKCkKZmlnCgpgYGAKCmBgYHtyfQojIyMgUGxvdCBlaWdlbmZ1bmN0aW9ucwoKb3B0aW9ucyhyZXByLnBsb3Qud2lkdGg9NiwgcmVwci5wbG90LmhlaWdodD0zKQpudW1iZXIubGFicyA8LSBwYXN0ZSgiRWlnZW5mdW5jdGlvbiIsIDE6bmVpZykKbmFtZXMobnVtYmVyLmxhYnMpIDwtIGMoIjEiOm5laWcpCmVpZ2VuZnVuY3Rpb25fdGliYmxlICU+JQogIGdncGxvdChhZXModGltZSwgdmFsdWUpKSArCiAgZ2VvbV9saW5lKGFlcyhsaW5ldHlwZSA9IGJvdW5kKSkgKwogIGZhY2V0X2dyaWQodHlwZSB+IG51bWJlciwgbGFiZWxsZXIgPSBsYWJlbGxlcihudW1iZXIgPSBudW1iZXIubGFicykpICsKICBzY2FsZV9saW5ldHlwZV9tYW51YWwodmFsdWVzPWMoImRhc2hlZCIsICJzb2xpZCIsICJkYXNoZWQiKSkrCiAgdGhlbWVfYncoKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiLCBzdHJpcC5iYWNrZ3JvdW5kID0gZWxlbWVudF9yZWN0KAogICAgICAgIGNvbG9yPSJibGFjayIsIGZpbGw9IiNGRkZGRkYiKSwgCiAgICAgICAgYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSA5MCwgdmp1c3QgPSAwLjUsIGhqdXN0PTEpKQogIApgYGAKCmBgYHtyfQojIyMgUGVyY2VudCB2YXJpYWJpbGl0eSBleHBsYWluZWQgYnkgZmlyc3QgZmV3IGVpZ2VuZnVuY3Rpb25zIGluIGZ1bmN0aW9uYWwgZGlyZWN0aW9uCjEwMCAqIHJvd01lYW5zKE1DTUNfZWlnZW4kZWlndmFsRnVuYylbcDE6KHAxLW5laWcrMSldCmBgYAoKYGBge3J9CiMjIyBQZXJjZW50IHZhcmlhYmlsaXR5IGV4cGxhaW5lZCBieSBmaXJzdCBmZXcgZWlnZW5mdW5jdGlvbnMgaW4gbG9uZ2l0dWRpbmFsIGRpcmVjdGlvbgoxMDAgKiByb3dNZWFucyhNQ01DX2VpZ2VuJGVpZ3ZhbExvbmcpW3AyOihwMi1uZWlnKzEpXQpgYGA=